# NOTE: you should use cleaned up data eventually here

#### Load Data
setwd("C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Data")  # Office computer directory
data <- get(load("ca_enrollment_data_AUG022019")) # individual-level data
households <- get(load("ca_household_data_AUG022019")) # household-level data
plan_data <- read.csv("ca_plan_data2.csv") # Covered California plan data
#network_data <- read.csv("ca_network_data2.csv",row.names=1) # CA network data
outside_logit <- get(load("sipp_logit"))
zip3_choices <- read.csv("zip3_choices2.csv",row.names = 1) # choice set by 3 digit zip and rating area
product_definitions <- read.csv("product_definitions.csv",row.names = 1) # definitions of column names in zip3_choices
library(plyr)

set.seed(1)

#### Determine whether household left market or just became uninsured

	# Add variables
		
		# Income brackets
		households$FPL_bracket <- "138orless"
		households[households$FPL > 1.38 & households$FPL <= 2.50 & !is.na(households$FPL),"FPL_bracket"] <- "138to250"
		households[households$FPL > 2.50 & households$FPL <= 4  & !is.na(households$FPL),"FPL_bracket"] <- "250to400"
		households[households$FPL > 4  & !is.na(households$FPL),"FPL_bracket"] <- "400ormore"
			
		# Age Brackets
		households$perc_18to34 <- households$perc_18to25 + households$perc_26to34
		households$perc_35to54 <- households$perc_35to44 + households$perc_45to54
		
		# Determine in/out of market
		households$out_of_market <- 0
		out_of_market_predictions <- as.numeric(predict(outside_logit,newdata=households[is.na(households$plan_number_nocsr),],type="response") >=
			runif(nrow(households[is.na(households$plan_number_nocsr),])))
		households[is.na(households$plan_number_nocsr),"out_of_market"] <- out_of_market_predictions
		
		# Remove records from market
		#households <- households[households$out_of_market == 0,]
		#data <- data[data$household_year %in% rownames(households),]

		# Record out of market in data object
		data$out_of_market <- households[data$household_year,"out_of_market"]

		# Get previous plan offered
		data$previous_plan_offered <- households[data$household_year,"previous_plan_offered"]

		# Delete household object
		rm(households)
		gc()
		
# NOTE: I am reordering plan data because I defined a plan_id number earlier based on a previous ordering
plan_data$Plan_ID_Order <- as.character(plan_data$Plan_ID_Order)
rownames(plan_data) <- paste(plan_data$Plan_Name,plan_data$ENROLLMENT_YEAR,
	plan_data$region,sep="")
plan_data <- plan_data[plan_data$Plan_ID_Order,]
		
		
plan_data$Issuer_Name <- as.character(plan_data$Issuer_Name)
plan_data[plan_data$Issuer_Name == "Anthem Blue Cross","Issuer_Name"] <- "Anthem"
plan_data[plan_data$Issuer_Name == "Blue Shield","Issuer_Name"] <- "Blue_Shield"
plan_data[plan_data$Issuer_Name == "Chinese Community","Issuer_Name"] <- "Chinese_Community"
plan_data[plan_data$Issuer_Name == "Contra Costa Health Plan","Issuer_Name"] <- "Contra_Costa"
plan_data[plan_data$Issuer_Name == "Health Net","Issuer_Name"] <- "Health_Net"
plan_data[plan_data$Issuer_Name == "LA Care","Issuer_Name"] <- "LA_Care"
plan_data[plan_data$Issuer_Name == "UnitedHealthcare","Issuer_Name"] <- "United"
plan_data[plan_data$Issuer_Name == "Western Health","Issuer_Name"] <- "Western"
plan_data[plan_data$Issuer_Name == "Sharp","Issuer_Name"] <- "SHARP"
	
# Create Age Group Variable

	data$age_group <- NA
	data[data$age < 18 & !is.na(data$age),"age_group"] <- "0to17"
	data[data$age >= 18 & data$age < 26 & !is.na(data$age),"age_group"] <- "18to25"
	data[data$age >= 26 & data$age < 35 & !is.na(data$age),"age_group"] <- "26to34"
	data[data$age >= 35 & data$age < 45 & !is.na(data$age),"age_group"] <- "35to44"
	data[data$age >= 45 & data$age < 55 & !is.na(data$age),"age_group"] <- "45to54"
	data[data$age >= 55 & data$age < 65 & !is.na(data$age),"age_group"] <- "55to64"
	data[data$age >= 65 & data$age < 120 & !is.na(data$age),"age_group"] <- "65plus"
	
# Metal

	plan_data$metal <- as.character(plan_data$metal_level)
	plan_data[plan_data$metal %in% c("Silver - Enhanced 73",
		"Silver - Enhanced 87","Silver - Enhanced 94"),"metal"] <- "Silver"
	
# Insurer-Metal

	plan_data$metal_short <- "CB"
	plan_data[plan_data$metal == "Silver","metal_short"] <- "SGP"
	#plan_data[plan_data$metal %in% c("Gold","Platinum"),"metal_short"] <- "ZGold_Plat"
	plan_data$insurer_metal_short <- paste(plan_data$metal_short,plan_data$Issuer_Name,sep="_")
	
# Market-year

	plan_data$year_market <- paste(plan_data$ENROLLMENT_YEAR,plan_data$region,sep="_")
	
# Insurer-Network Type

	plan_data$HMO <- "PPO"
	plan_data[plan_data$PLAN_NETWORK_TYPE == "HMO","HMO"] <- "HMO"

	plan_data$insurer_net <- paste(plan_data$Issuer_Name,plan_data$HMO,sep="_")
	plan_data$insurer_net_metal <- paste(plan_data$metal_short,plan_data$insurer_net,sep="_")
	#plan_data$insurer_net_metal <- paste(plan_data$Issuer_Name,plan_data$PLAN_NETWORK_TYPE,plan_data$metal,sep="_")
	
# Rank plans by premiums
	
# End Week

	# Important: a number of records have missing values for end week
	# I think it's likely that the end week is the end of the year
	
	data[is.na(data$end_week),"end_week"] <- 52
	gc()
	
# Functions

	check_whether_plan_offered <- function(i,zip3_choices,previous_plans) {
		
		region_year_plans <- which(plan_data$ENROLLMENT_YEAR == zip3_choices[i,"Year"] &
			plan_data$region == zip3_choices[i,"Region"])
		
		# Get the broad product categories available (insurer/network type/etc. - no metal)
		available_products <- zipchoices[i,]
		available_products <- names(available_products)[!is.na(available_products)]
		available_insurers <- unique(product_definitions[available_products,"insurer"])
			
			# Now get the specific plans available
			
			available_plans <- c()
			if ("Anthem" %in% available_insurers) {
				if("Anthem_HMO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Anthem" & plan_data$PLAN_NETWORK_TYPE == "HMO" &
							plan_data$MSP == 0)))
				}
				if("Anthem_EPO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Anthem" & plan_data$PLAN_NETWORK_TYPE == "EPO" &
							plan_data$MSP == 0)))
				}
				if("Anthem_PPO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Anthem" & plan_data$PLAN_NETWORK_TYPE == "PPO" &
							plan_data$MSP == 0)))
				}
				if("Anthem_EPO_MSP" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Anthem" & plan_data$PLAN_NETWORK_TYPE == "EPO" &
							plan_data$MSP == 1)))
				}
				if("Anthem_PPO_MSP" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Anthem" & plan_data$PLAN_NETWORK_TYPE == "PPO" &
							plan_data$MSP == 1)))
				}
			} 
			if ("Blue_Shield" %in% available_insurers) {
				# could be HMO, EPO, PPO
				if("Blue_Shield_HMO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Blue_Shield" & plan_data$PLAN_NETWORK_TYPE == "HMO" )))
				}
				if("Blue_Shield_EPO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Blue_Shield" & plan_data$PLAN_NETWORK_TYPE == "EPO" )))
				}
				if("Blue_Shield_PPO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Blue_Shield" & plan_data$PLAN_NETWORK_TYPE == "PPO" )))
				}
			} 
			if ("Health_Net" %in% available_insurers) {
				# could be HMO, HSP, EPO, PPO
				if("Health_Net_HMO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Health_Net" & plan_data$PLAN_NETWORK_TYPE == "HMO")))
				}
				if("Health_Net_HSP" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Health_Net" & plan_data$PLAN_NETWORK_TYPE == "HSP")))
				}
				if("Health_Net_EPO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Health_Net" & plan_data$PLAN_NETWORK_TYPE == "EPO")))
				}
				if("Health_Net_PPO" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "Health_Net" & plan_data$PLAN_NETWORK_TYPE == "PPO")))
				}
			}
			if ("SHARP" %in% available_insurers) {
				# could be network 1 or 2
				if("Sharp1" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "SHARP" & plan_data$Network_num == 1 & 
							!is.na(plan_data$Network_num))))
				}
				if("Sharp2" %in% available_products) {
					available_plans <- c(available_plans,
						intersect(region_year_plans,
							which(plan_data$Issuer_Name == "SHARP" & plan_data$Network_num == 2 & 
							!is.na(plan_data$Network_num))))
				}
			}
			if (any(single_product_insurers %in% available_insurers)) {
				available_plans <- c(available_plans,
					intersect(region_year_plans,
						which(plan_data$Issuer_Name %in% intersect(available_insurers,single_product_insurers))))
			}	
			
		# Now see if previous plan in set of available plans
		available_plans <- plan_data[available_plans,"Plan_Name_NOCSR"]
		plan_offered <- previous_plans %in% available_plans
		return(plan_offered)
	}

	get_plans <- function(transdata,year1,year2) {
	
		participants_year1 <- which(transdata$year == year1)
		choice_year1 <- transdata[participants_year1,"plan_id"]
		
		id_years_to_check <- paste(transdata[participants_year1,"individual_id"],year2,sep="_")
		returnees <- id_years_to_check[id_years_to_check %in% rownames(transdata)]
		choice_year2 <- rep(NA,length(id_years_to_check))
		names(choice_year2) <- id_years_to_check
		choice_year2[returnees] <- transdata[returnees,"plan_id"]
		choices <- data.frame(choice_year1,choice_year2)
		choices[,1] <- as.character(choices[,1])
		choices[,2] <- as.character(choices[,2])
		
		# Add Demographic variables
		choices$OEP <- transdata[participants_year1,"OEP"]
		choices$age_group <- transdata[participants_year1,"age_group"]
		choices$region <- transdata[participants_year1,"region"]
		choices$zip_region_year <- paste(transdata[participants_year1,"zip_3"],
			transdata[participants_year1,"region"],year2,sep="_")
		choices$flagged <- transdata[participants_year1,"flagged"]
		choices$end_week <- transdata[participants_year1,"end_week"]
		
		choices$end_week_2017 <- NA
		choices[returnees,"end_week_2017"] <- transdata[returnees,"end_week"]
		
		return(choices)
	}
	
	get_new_consumers <- function(transdata,year1,year2) {
	
		participants_year2 <- which(transdata$year == year2)
		choice_year2 <- transdata[participants_year2,"plan_id"]
		
		id_years_to_check <- paste(transdata[participants_year2,"individual_id"],year1,sep="_")
		returnees <- id_years_to_check[id_years_to_check %in% rownames(transdata)]
		choice_year1 <- rep(NA,length(id_years_to_check))
		names(choice_year1) <- id_years_to_check
		choice_year1[returnees] <- transdata[returnees,"plan_id"]
		choices <- data.frame(choice_year1,choice_year2)
		choices[,1] <- as.character(choices[,1])
		choices[,2] <- as.character(choices[,2])
		
		# Add Demographic variables
		choices$OEP <- transdata[participants_year2,"OEP"]
		choices$age_group <- transdata[participants_year2,"age_group"]
		choices$region <- transdata[participants_year2,"region"]
		choices$zip_region_year <- paste(transdata[participants_year2,"zip_3"],
			transdata[participants_year2,"region"],year2,sep="_")
		choices$flagged <- transdata[participants_year2,"flagged"]
		choices$end_week <- transdata[participants_year2,"end_week"]
		
		choices$end_week_2017 <- NA
		choices[returnees,"end_week_2017"] <- transdata[returnees,"end_week"]
	
		return(choices)
	}
	
	# NOTE: should you remove flagged people from this?
	
	tabulate_output <- function(planids,plans,insurers) {
	
		departed_before <- which(is.na(plans[,"choice_year2"]) & 
			planids[,"end_week"] < 52 & !is.na(planids[,"end_week"]))
		departed_after <- which(is.na(plans[,"choice_year2"]) & 
			planids[,"end_week"] >= 52 & !is.na(planids[,"end_week"]))
		departed <- which(is.na(plans[,"choice_year2"]))
		not_offered <- setdiff(which(!planids$plan_offered),departed)
		not_offered_same <- length(intersect(not_offered,
			which(plans[,"choice_year1"] == plans[,"choice_year2"])))
		not_offered_diff <- length(intersect(not_offered,
			which(plans[,"choice_year1"] != plans[,"choice_year2"])))
	
		# NOTE: 112 records supposedly remained on the same plan between 2016 and 2017, but their plan 
		# was no longer offered in 2017 - I'm counting them in not offered (too small to worry about)
		
		same <- length(setdiff(which(plans[,"choice_year1"] == plans[,"choice_year2"]),not_offered))
		diff_plan_same_ins <- length(setdiff(which(plans[,"choice_year1"] != plans[,"choice_year2"] &
			!is.na(plans[,"choice_year2"]) & 
			insurers[,"choice_year1"] == insurers[,"choice_year2"]),not_offered))
		diff_plan_diff_ins <- length(setdiff(which(plans[,"choice_year1"] != plans[,"choice_year2"] &
			!is.na(plans[,"choice_year2"]) & 
			insurers[,"choice_year1"] != insurers[,"choice_year2"]),not_offered))
		not_offered <- length(not_offered)
		
		departed_before <- length(departed_before)
		departed_after <- length(departed_after)
		departed_offered <- length(intersect(departed,which(planids$plan_offered)))
		departed_not_offered <- length(intersect(departed,which(!planids$plan_offered)))
		departed <- length(departed)
		
		output <- c(same,not_offered,not_offered_same,not_offered_diff,
			diff_plan_same_ins,diff_plan_diff_ins,departed,
			departed_before,departed_after,departed_offered,
			departed_not_offered)
	
		return(c(output,sum(output,na.rm=TRUE)))
	}
	
	
# Recreate dynamic choices

	individual_names <- unique(data$individual_id)
	years <- c(2014:2018)
	dynamic_choices <- matrix(NA,length(individual_names),length(years), dimnames=list(individual_names,years))
	prevoffer_choices <- matrix(NA,length(individual_names),length(years), dimnames=list(individual_names,years))
	markets_choices <- matrix(NA,length(individual_names),length(years), dimnames=list(individual_names,years))
			
	for(t in years) {
		names_t <- as.character(data[data$year == t & !data$flagged,"individual_id"])				
		plan_numbers_year <- as.numeric(data[data$year == t & !data$flagged,"plan_id"])
		markets_year <- as.numeric(data[data$year == t & !data$flagged,"rating_area"])
		
		dynamic_choices[names_t,as.character(t)] <- plan_numbers_year
		markets_choices[names_t,as.character(t)] <- markets_year
		
		if(t==2014) {
			prevoffer_choices[names_t,as.character(t)] <- plan_numbers_year
		} else {
			prevoffer_names_t <- as.character(data[data$year == t & !data$flagged & 
				data$previous_plan_offered == 1 & !is.na(data$previous_plan_offered),"individual_id"])				
			prevoffer_plan_numbers_year <- as.numeric(data[data$year == t & !data$flagged & 
				data$previous_plan_offered == 1 & !is.na(data$previous_plan_offered),"plan_id"])
			prevoffer_choices[prevoffer_names_t,as.character(t)] <- prevoffer_plan_numbers_year
		}
		

	}		
	gc()
	
	colnames(markets_choices) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	plan_data$Plan_Name2_NOCSR <- as.character(plan_data$Plan_Name2_NOCSR)
	plan_trans <- cbind(plan_data[dynamic_choices[,"2014"],"Plan_Name2_NOCSR"],
		plan_data[dynamic_choices[,"2015"],"Plan_Name2_NOCSR"],
		plan_data[dynamic_choices[,"2016"],"Plan_Name2_NOCSR"],
		plan_data[dynamic_choices[,"2017"],"Plan_Name2_NOCSR"],
		plan_data[dynamic_choices[,"2018"],"Plan_Name2_NOCSR"])
	rownames(plan_trans) <- rownames(dynamic_choices)
	colnames(plan_trans) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
		# Tenure in Exchange
		tenure <- matrix(NA,nrow(plan_trans),ncol(plan_trans),dimnames=list(rownames(plan_trans),colnames(plan_trans)))
		
		tenure[,"Y2014"] <- as.numeric(!is.na(plan_trans[,"Y2014"]))
		tenure[,"Y2015"] <- as.numeric(!is.na(plan_trans[,"Y2015"]))
		tenure[,"Y2016"] <- as.numeric(!is.na(plan_trans[,"Y2016"]))
		tenure[,"Y2017"] <- as.numeric(!is.na(plan_trans[,"Y2017"]))
		tenure[,"Y2018"] <- as.numeric(!is.na(plan_trans[,"Y2018"]))
		
		tenure[tenure[,"Y2015"] > 0,"Y2015"] <- tenure[tenure[,"Y2015"] > 0,"Y2014"] + 1
		tenure[tenure[,"Y2016"] > 0,"Y2016"] <- tenure[tenure[,"Y2016"] > 0,"Y2015"] + 1
		tenure[tenure[,"Y2017"] > 0,"Y2017"] <- tenure[tenure[,"Y2017"] > 0,"Y2016"] + 1
		tenure[tenure[,"Y2018"] > 0,"Y2018"] <- tenure[tenure[,"Y2018"] > 0,"Y2017"] + 1
		
		plan_tenure <- matrix(NA,nrow(plan_trans),ncol(plan_trans),dimnames=list(rownames(plan_trans),colnames(plan_trans)))
		
		plan_tenure[,"Y2014"] <- as.numeric(!is.na(plan_trans[,"Y2014"]))
		plan_tenure[,"Y2015"] <- as.numeric(!is.na(plan_trans[,"Y2015"]))
		plan_tenure[,"Y2016"] <- as.numeric(!is.na(plan_trans[,"Y2016"]))
		plan_tenure[,"Y2017"] <- as.numeric(!is.na(plan_trans[,"Y2017"]))
		plan_tenure[,"Y2018"] <- as.numeric(!is.na(plan_trans[,"Y2018"]))

		
		plan_tenure[plan_tenure[,"Y2014"] > 0 & plan_tenure[,"Y2015"] > 0 & (plan_trans[,"Y2014"] == plan_trans[,"Y2015"]),"Y2015"] <- 
			plan_tenure[plan_tenure[,"Y2014"] > 0 & plan_tenure[,"Y2015"] > 0 & (plan_trans[,"Y2014"] == plan_trans[,"Y2015"]),"Y2014"] + 1
		plan_tenure[plan_tenure[,"Y2015"] > 0 & plan_tenure[,"Y2016"] > 0 & (plan_trans[,"Y2015"] == plan_trans[,"Y2016"]),"Y2016"] <- 
			plan_tenure[plan_tenure[,"Y2015"] > 0 & plan_tenure[,"Y2016"] > 0 & (plan_trans[,"Y2015"] == plan_trans[,"Y2016"]),"Y2015"] + 1
		plan_tenure[plan_tenure[,"Y2016"] > 0 & plan_tenure[,"Y2017"] > 0 & (plan_trans[,"Y2016"] == plan_trans[,"Y2017"]),"Y2017"] <- 
			plan_tenure[plan_tenure[,"Y2016"] > 0 & plan_tenure[,"Y2017"] > 0 & (plan_trans[,"Y2016"] == plan_trans[,"Y2017"]),"Y2016"] + 1
		plan_tenure[plan_tenure[,"Y2017"] > 0 & plan_tenure[,"Y2018"] > 0 & (plan_trans[,"Y2017"] == plan_trans[,"Y2018"]),"Y2018"] <- 
			plan_tenure[plan_tenure[,"Y2017"] > 0 & plan_tenure[,"Y2018"] > 0 & (plan_trans[,"Y2017"] == plan_trans[,"Y2018"]),"Y2017"] + 1
		
		plan_tenure_2016 <- table(plan_tenure[plan_tenure[,"Y2016"] > 0,"Y2016"])/sum(table(plan_tenure[plan_tenure[,"Y2016"] > 0,"Y2016"]))
		plan_tenure_2017 <- table(plan_tenure[plan_tenure[,"Y2017"] > 0,"Y2017"])/sum(table(plan_tenure[plan_tenure[,"Y2017"] > 0,"Y2017"]))
		plan_tenure_2018 <- table(plan_tenure[plan_tenure[,"Y2018"] > 0,"Y2018"])/sum(table(plan_tenure[plan_tenure[,"Y2018"] > 0,"Y2018"]))

		plan_tenure_output <- cbind(c(plan_tenure_2016,NA,NA),c(plan_tenure_2017,NA),plan_tenure_2018)
		write.csv(plan_tenure_output,"plan_tenure_output.csv")
		
	prevoffer_plan_trans <- cbind(plan_data[prevoffer_choices[,"2014"],"Plan_Name2_NOCSR"],
		plan_data[prevoffer_choices[,"2015"],"Plan_Name2_NOCSR"],
		plan_data[prevoffer_choices[,"2016"],"Plan_Name2_NOCSR"],
		plan_data[prevoffer_choices[,"2017"],"Plan_Name2_NOCSR"],
		plan_data[prevoffer_choices[,"2018"],"Plan_Name2_NOCSR"])
	rownames(prevoffer_plan_trans) <- rownames(prevoffer_choices)
	colnames(prevoffer_plan_trans) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	
	insurers_trans <- cbind(plan_data[dynamic_choices[,"2014"],"Issuer_Name"],
		plan_data[dynamic_choices[,"2015"],"Issuer_Name"],
		plan_data[dynamic_choices[,"2016"],"Issuer_Name"],
		plan_data[dynamic_choices[,"2017"],"Issuer_Name"],
		plan_data[dynamic_choices[,"2018"],"Issuer_Name"])
	rownames(insurers_trans) <- rownames(dynamic_choices)
	colnames(insurers_trans) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	insurers_net_trans <- cbind(plan_data[dynamic_choices[,"2014"],"insurer_net"],
		plan_data[dynamic_choices[,"2015"],"insurer_net"],
		plan_data[dynamic_choices[,"2016"],"insurer_net"],
		plan_data[dynamic_choices[,"2017"],"insurer_net"],
		plan_data[dynamic_choices[,"2018"],"insurer_net"])
	rownames(insurers_net_trans) <- rownames(dynamic_choices)
	colnames(insurers_net_trans) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	metals_trans <- cbind(plan_data[dynamic_choices[,"2014"],"metal"],
		plan_data[dynamic_choices[,"2015"],"metal"],
		plan_data[dynamic_choices[,"2016"],"metal"],
		plan_data[dynamic_choices[,"2017"],"metal"],
		plan_data[dynamic_choices[,"2018"],"metal"])
	rownames(metals_trans) <- rownames(dynamic_choices)
	colnames(metals_trans) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	insurer_metal_trans <- cbind(plan_data[dynamic_choices[,"2014"],"insurer_metal_short"],
		plan_data[dynamic_choices[,"2015"],"insurer_metal_short"],
		plan_data[dynamic_choices[,"2016"],"insurer_metal_short"],
		plan_data[dynamic_choices[,"2017"],"insurer_metal_short"],
		plan_data[dynamic_choices[,"2018"],"insurer_metal_short"])
	rownames(insurer_metal_trans) <- rownames(dynamic_choices)
	colnames(insurer_metal_trans) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	insurer_net_metal_trans <- cbind(plan_data[dynamic_choices[,"2014"],"insurer_net_metal"],
		plan_data[dynamic_choices[,"2015"],"insurer_net_metal"],
		plan_data[dynamic_choices[,"2016"],"insurer_net_metal"],
		plan_data[dynamic_choices[,"2017"],"insurer_net_metal"],
		plan_data[dynamic_choices[,"2018"],"insurer_net_metal"])
	rownames(insurer_net_metal_trans) <- rownames(dynamic_choices)
	colnames(insurer_net_metal_trans) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	year_labels <- c(2014:2018)
	names(year_labels) <- c("Y2014","Y2015","Y2016","Y2017","Y2018")
	
	for(t in year_labels) {
		year_label_t <- names(which(year_labels == t))
		out_indices <- which(is.na(insurers_trans[,year_label_t]))
		out_individuals <- paste(individual_names[out_indices],t,sep="_")
		out_individual_indices <- which(out_individuals %in% rownames(data))
		out_individuals <- out_individuals[out_individual_indices]
		
		went_to_outside_option <- out_individual_indices[which(data[out_individuals,"out_of_market"]  == 0)]
		out_of_market_class <- rep("Exited",length(out_indices))
		out_of_market_class[went_to_outside_option] <- "Outside"
		
		plan_trans[out_indices,year_label_t] <- out_of_market_class
		insurers_trans[out_indices,year_label_t] <- out_of_market_class
		insurers_net_trans[out_indices,year_label_t] <- out_of_market_class
		metals_trans[out_indices,year_label_t] <- out_of_market_class
		insurer_metal_trans[out_indices,year_label_t] <- out_of_market_class
		insurer_net_metal_trans[out_indices,year_label_t] <- out_of_market_class
	}
	
#### Transition Tables

	# Overall plans
	trans_plan_names <- sort(unique(plan_data$Plan_Name2_NOCSR))
	transitions <- matrix(0,length(trans_plan_names),length(trans_plan_names),dimnames=list(trans_plan_names,trans_plan_names))
	#prevoffer_transitions <- matrix(0,length(trans_plan_names),length(trans_plan_names),dimnames=list(trans_plan_names,trans_plan_names))
	
	transitions1415 <- xtabs(formula = ~ Y2014 + Y2015,
		data=plan_trans[!(is.na(plan_trans[,"Y2014"]) & is.na(plan_trans[,"Y2015"])) & 
			!(is.na(markets_choices[,"Y2014"]) & is.na(markets_choices[,"Y2015"])),])
	transitions1415 <- transitions1415[setdiff(rownames(transitions1415),c("Outside","Exited")),
		setdiff(colnames(transitions1415),c("Outside","Exited"))]
	transitions[rownames(transitions1415),colnames(transitions1415)] <- 
		transitions[rownames(transitions1415),colnames(transitions1415)] + transitions1415
	
	transitions1516 <- xtabs(formula = ~ Y2015 + Y2016,
		data=plan_trans[!(is.na(plan_trans[,"Y2015"]) & is.na(plan_trans[,"Y2016"])) & 
			!(is.na(markets_choices[,"Y2015"]) & is.na(markets_choices[,"Y2016"])),])	
	transitions1516 <- transitions1516[setdiff(rownames(transitions1516),c("Outside","Exited")),
		setdiff(colnames(transitions1516),c("Outside","Exited"))]
	transitions[rownames(transitions1516),colnames(transitions1516)] <- 
		transitions[rownames(transitions1516),colnames(transitions1516)] + transitions1516
	
	transitions1617 <- xtabs(formula = ~ Y2016 + Y2017,
		data=plan_trans[!(is.na(plan_trans[,"Y2016"]) & is.na(plan_trans[,"Y2017"])) & 
			!(is.na(markets_choices[,"Y2016"]) & is.na(markets_choices[,"Y2017"])),])	
	transitions1617 <- transitions1617[setdiff(rownames(transitions1617),c("Outside","Exited")),
		setdiff(colnames(transitions1617),c("Outside","Exited"))]
	transitions[rownames(transitions1617),colnames(transitions1617)] <- 
		transitions[rownames(transitions1617),colnames(transitions1617)] + transitions1617
	
	transitions1718 <- xtabs(formula = ~ Y2017 + Y2018,
		data=plan_trans[!(is.na(plan_trans[,"Y2017"]) & is.na(plan_trans[,"Y2018"])) & 
			!(is.na(markets_choices[,"Y2017"]) & is.na(markets_choices[,"Y2018"])),])	
	transitions1718 <- transitions1718[setdiff(rownames(transitions1718),c("Outside","Exited")),
		setdiff(colnames(transitions1718),c("Outside","Exited"))]
	transitions[rownames(transitions1718),colnames(transitions1718)] <- 
		transitions[rownames(transitions1718),colnames(transitions1718)] + transitions1718
	
	transitions <- transitions/4
	
	sum(diag(transitions))/sum(transitions)
	
	# Overall - Metals

	metal_order <- c("Minimum Coverage","Bronze","Silver","Gold","Platinum","Outside","Exited")
	transitions1415 <- xtabs(formula = ~ Y2014 + Y2015,
		data=metals_trans[!(is.na(metals_trans[,"Y2014"]) & is.na(metals_trans[,"Y2015"])) & 
			!(is.na(markets_choices[,"Y2014"]) & is.na(markets_choices[,"Y2015"])),])	
	transitions1516 <- xtabs(formula = ~ Y2015 + Y2016,
		data=metals_trans[!(is.na(metals_trans[,"Y2015"]) & is.na(metals_trans[,"Y2016"])) & 
			!(is.na(markets_choices[,"Y2015"]) & is.na(markets_choices[,"Y2016"])),])	
	transitions1617 <- xtabs(formula = ~ Y2016 + Y2017,
		data=metals_trans[!(is.na(metals_trans[,"Y2016"]) & is.na(metals_trans[,"Y2017"])) & 
			!(is.na(markets_choices[,"Y2016"]) & is.na(markets_choices[,"Y2017"])),])	
	transitions1718 <- xtabs(formula = ~ Y2017 + Y2018,
		data=metals_trans[!(is.na(metals_trans[,"Y2017"]) & is.na(metals_trans[,"Y2018"])) & 
			!(is.na(markets_choices[,"Y2017"]) & is.na(markets_choices[,"Y2018"])),])	
	transitions <- rbind(transitions1415[c("Minimum Coverage","Bronze","Silver","Gold","Platinum"),metal_order],0,0)	+ 
		transitions1516[metal_order,metal_order] + transitions1617[metal_order,metal_order] + 
		transitions1718[metal_order,metal_order]
	
	rownames(transitions) <- metal_order
	transitions_perc <- transitions/rowSums(transitions)
	
	write.csv(transitions,"output.csv")
	write.csv(transitions/rowSums(transitions),"output_perc.csv")
	
	# Overall - Insurers

	big_insurers <- c("Anthem","Blue_Shield","Health_Net","Kaiser")
	short_insurers_trans <- insurers_trans
	for(t in year_labels) {
		year_label_t <- names(which(year_labels == t))
		short_insurers_trans[!short_insurers_trans[,year_label_t] %in% 
			c(big_insurers,"Outside","Exited"),year_label_t] <- "Other_Insurer"
	}
	
	insurer_order <- c("Anthem","Blue_Shield","Health_Net","Kaiser","Other_Insurer","Outside","Exited")
	transitions1415 <- xtabs(formula = ~ Y2014 + Y2015,
		data=short_insurers_trans[!(is.na(short_insurers_trans[,"Y2014"]) & is.na(short_insurers_trans[,"Y2015"])) & 
			!(is.na(markets_choices[,"Y2014"]) & is.na(markets_choices[,"Y2015"])),])	
	transitions1516 <- xtabs(formula = ~ Y2015 + Y2016,
		data=short_insurers_trans[!(is.na(short_insurers_trans[,"Y2015"]) & is.na(short_insurers_trans[,"Y2016"])) & 
			!(is.na(markets_choices[,"Y2015"]) & is.na(markets_choices[,"Y2016"])),])	
	transitions1617 <- xtabs(formula = ~ Y2016 + Y2017,
		data=short_insurers_trans[!(is.na(short_insurers_trans[,"Y2016"]) & is.na(short_insurers_trans[,"Y2017"])) & 
			!(is.na(markets_choices[,"Y2016"]) & is.na(markets_choices[,"Y2017"])),])	
	transitions1718 <- xtabs(formula = ~ Y2017 + Y2018,
		data=short_insurers_trans[!(is.na(short_insurers_trans[,"Y2017"]) & is.na(short_insurers_trans[,"Y2018"])) & 
			!(is.na(markets_choices[,"Y2017"]) & is.na(markets_choices[,"Y2018"])),])	
	transitions <- rbind(transitions1415[c("Anthem","Blue_Shield","Health_Net","Kaiser","Other_Insurer"),insurer_order],0,0)	+ 
		transitions1516[insurer_order,insurer_order] + transitions1617[insurer_order,insurer_order] + 
		transitions1718[insurer_order,insurer_order]
	
	rownames(transitions) <- insurer_order
	transitions_perc <- transitions/rowSums(transitions)
	
	write.csv(transitions,"output.csv")
	write.csv(transitions/rowSums(transitions),"output_perc.csv")
	
